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Abstract 

The most difficult aspect of the realistic modeling of granular materials is how 
to capture the real shape of the particles. Here we present a method to simulate 
granular materials with complex-shaped particles. The particle shape is represented 
by the classical concept of a Minkowski sum, which permits the representation of 
complex shapes without the need to define the object as a composite of spherical or 
convex particles. A well defined interaction force between these bodies is derived. 
The algorithm for identification of neighbor particles reduces force calculations to 
O(N), where ./V is the number of particles. The algorithm is much more efficient, 
accurate and easier to implement than other models.. We investigate the existence 
of a statistical equilibrium in granular systems with circular non-spherical particles 
in the collisional. regime. We also investigate the limit state of dissipative granular 
materials using biaxial test simulations. The results agree with the classical assump- 
tion of the statistical mechanics for non-dissipative systems, and the critical state 
theory of soils mechanics for dissipative granular materials. 

Key words: Granular systems, Dynamics and kinematics of rigid bodies, 
Molecular dynamics methods 
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1 Introduction 



Rapid advances in computer simulations have led to many new developments 
in the modeling of particulate systems. These systems represent different real 
physical systems at different scales, such as the small scale of liquid crystals 
PPj, geological scales of snow and debris flow, and the astronomical scales of 
planetary rings or dynamics evolution of precursors of planets [2]. Although 
particle shape plays an important role, most theoretical and numerical de- 
velopments have been restricted to particles with spherical or circular shape. 
These simplification lead in some cases to unrealistic properties. In collisional 



non-dissipative systems, spherical (or circular) particles can not exchange an- 
gular momentum, so that the system cannot explore all the phase space dur- 
ing the evolution. In dissipative granular systems such as sand piles or fault 
gouges, disks of spheres tend to roll more easily than non spherical particles, 
leading to unrealistic angles of repose and bulk friction coefficients. 

Three different approach has been presented to model the real shape of par- 
ticulate materials. In the first approach the shape is represented as a scalar 
funcions. This model allows to represent especific shapes, such as ellipses [3], 
ellipsoids [4 J and superquadric [3] . The main drawback of those methods is that 
the calculations required in the contact force are much more expensive than 
in spherical (or circular) particles. In the second approach the non-spherical 
particle is represented ad aggregates or clumps of disks and spheres bonded 
together [5117] . In this approach crushing and fracture of aggregates can be eas- 
ily modeled. The disadvantage is that this method requires a larger number 
of particles. 

The third approach is represent the complex shape using polygons in 2D, 
[51I51TTU] or polyhedrons in 3D [TlfTS] The most difficult aspect for the simula- 
tions of these objects is the handling of contact interactions. Nowadays, the 
interaction is resolved by decomposing them in convex pieces, and applying 
penalty methods, impulse-based methods or dynamic constraints in the inter- 
action between these pieces. Impulse-based methods (also called even-driven 
methods) allow real-time simulations, but they cannot handle permanent or 
lasting contacts [TT]. On the other hand, constraint methods (or contact dy- 
namics methods) can handle resting contacts with infinite stiffness, but simula- 
tions are computationally expensive and lead in some cases to indeterminacy 
in the solution of contact forces . This indeterminacy is removed by using 
penalty methods, where the bodies are allowed to interpenetrate each other 
and the force is calculated in terms of their overlap. However, until very re- 
cently the determination of such contact force has been heuristic and lacks 
physical correctness, because the interactions do not comply with energy con- 
servation [13] 

We propose a solution to this problem in 2D based in the mathematical concept 
of spheropolygons. These objects are generated from the Minkowski sum of 
a polygon with a disk, which is nothing more than the object resulting from 
sweeping a disk around the polygon. The 3D counterpart of these objects 
are the spheropolyhedrons, which were recently introduced by Pournin and 
Liebling for the simulations of granular media. This simple concept can be 
used to generate very complex shapes, including non-convex bodies, without 
the need to decompose them into spherical or convex parts. Here we introduce 
an efficient method to calculate dissipative and non-dissipative interaction 
between spheropolygons. We probe also that the model complies with the 
statistical mechanical principles and the physical laws of a conservative system. 
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Since the code can simulate particles a wide range of shapes, it can be used to 
investigate the effect of the aspect ratio, angularity and non-convexity of the 
particles on granular flow. The relevance of this investigation is demonstrated 
in the stress-strain response in biaxial test bellow, where we show that the 
system reach the critical state of soil mechanics, where the particle shape 
strongly affect the material properties of the granular media. 



2 Model 



Systems with different particle shapes are modeled using the concept of the 
Minkowski sum of a polygon (or polyline) with a disk. This mathematical 
operation is explained in Sub-Section 12. 1[ Sub-section 12.21 leads with the nu- 
merical calculation of mass properties . The interaction between the particles 
is obtained from the individual interaction between each vertex of one polygon 
and each edge of another, as explained in Sub-Section 12.31 In Sub-Section 12.41 
we present a method to redocue the number of floating point operations used 
to calculate interaction forces is drastically reduced by using a neighbor list 
and a contact list for each pair of neighbor particles. In Sub-Section 12.51 we 
present the algorith used in each time step of the dynamic simulation. 



2.1 Minkowski sum 



Given two sets of points P and Q in an Euclidean space, their Minkowski sum 
is given by P + Q = {x + y\ x&P, y G Q}. This operation is geometrically 
equivalent to the sweeping of one set around the profile of the other without 
changing the relative orientation. A special case is the sum of a polygon with a 
disk, which is defined here as spheropolygon. Other examples of a Minkowski 
sum are the spherocyllinder (sphere + line segment) [14|, the spherosimplex 
(sphere-l- simplex) [15] and the spheropolyhedron (sphere+polyhedron) |16j . 
which are used in simulations of particulate systems. 

The main advantage of the spheropolygons is that they allow us to represent 
any shape in 2D, from rounded to angular particles, and from convex to non- 
convex shapes. As we will see in Sub-Section 12.31 The Minkowsky sum does 
not need to be explicitely calculated to determine the particle interaction. 
The calculation of the mass properties, however, are calculated numerically, 
but this does not affect much the simulation time because the calculations are 
done only a the beginning of the simulations. 
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Fig. 1. Minkowski sum of a polygon with a disk. 

2.2 Mass properties 

Before calculating the mass, center of mass and moment of inertica of the 
spheropolygons we need to introduce some useful concepts. Given a spheropoly- 
gon SP = P+S, the polygon P will be called the polygon-base; and the radius 
r of the disk S sphero-radius. The distance d(X, P) from a test point X to 
the polygon-base is defined as follows: If the point is outside the polygon, it 
is given by the minimum distance between the point and the edges of the 
polygon; If the point is inside the polygon, we let d(X, P) = 0. Finally, the 
point X is inside of the spheropolygon when it satisfies d(X, P) < r. 

The point-inside-spheropolygon test combined with a basic Monte Carlo method 
is used to evaluate the integral expressions for mass, center of mass and the 
moment of inertia. The numerical integration is performed by taking a quasi- 
random set of points Xi uniformly distributed in a rectangular box containing 
the object. Then the integral over the area enclosed by the spheropolygon of 
any function j{X) is calculated as: 

r A Np 

M= J f(X)da*^Y.x(X)f(X)- (1) 



Ab ox is the area of the rectangular box, Xi is a quasi-random point inside Ab ox , 
and N p = 1.6 x 10 4 is the number of points. x{X) is the characteristic func- 
tion, which returns one if X is inside the spheropolygon and zero otherwise. 
Replacing f(X) by a, aX <r||X|| 2 results in M = m,mf,I + m||r|| 2 respec- 
tively, where a is the density, and m, r and / are the mass, center of mass and 
moment of inertia. 
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2.3 Interaction force 



To solve the interaction between spheropolygons we consider all vertex-edge 
distances between the polygons base, we consider two spheropolygons SPi 
and SPj with their respective polygons base P; and Pj and sphero-radii and 
Tj. Each polygon is defined by the set of vertices {Vj} and edges {Ej}. The 
overlapping length between each pair of vertex-edge (Vj, Ej) is defined as 

5(V,E) = {n + rj -d(V,E)), (2) 



where d(X, E) = \ \Y — X\ \ is the Euclidean distance from the vertex V to the 
segment E. Here X is the position of the vertex V and Y is its closest point on 
the edge E. The ramp function (x) returns x if x > and zero otherwise. The 
overlapping length in Eq. (jSJ) is equivalent to the interpenetration between the 
disks of radii rj and rj centered on X and Y. 

The force Fj acting on particle % by the particle j is defined by: 

4 = -Fa = E F(Yi, Ej) + E ^ ^<)> (3) 

Vi& VjE t 



where F(V, E) represent the force between the vertex V and the edge E. 
if the vertex-edge pair do not overlap, F(V, E) = 0. Different of vertex-edge 
forces can be included in the modelL linear dashpots, non-linear Hertzian laws, 
dissipative viscous forces proportional to the relative normal and tangential 
velocities, sliding friction, etc The force of Eq. [3] is applied to each particle in 
the middle point of the overlap region between the vertex and the edge: 

R{V, E)=X+(r l + U{V, E)) *~ Y (4) 

^ \\Y — -A 



so that the resulting torque on particle i given by j is 

nj = YsViEj (R{V h Ej) - rj) x F(E h Vj) 

^ (5) 
+ Y.v jEi (R(Vj, Ei) - fj x F{Ej, Vi), 



where is the center of mass of particle i. 

The evolution of f\ and the orientation tpi of the particle is governed by the 
equations of motion: 

rriifi = Fji ~ m iQy, Ii<Pi = E T ji- ( 6 ) 
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Here mi and Ij are the mass and moment of inertia of the particle. The sum 
is over all particles interacting with this particle; g is the gravity; and y is the 
unit vector along the vertical direction. 



2.4 Efficient calculation of forces 

The efficiency of the dynamics simulation is mainly determined by the method 
of contact detection. In a system of N particles, each one with M edges, the 
number of operations required to update the positions of the particle in each 
time step is in the order of O(NM), whereas the number of calculations for 
contact detection is 0(N 2 M 2 ). Simulations therefore become very slow when 
either the number of particles or the number of vertices is large. 

The first step to speed up the simulations is to execute the force calculation 
only over neighbor particles. With this aim we introduce the neighbor list, 
which is the collection of pair particles whose distance between them is less 
than 25. (The distance between two particles is defined as the minimum of all 
vertex-edge distances). The parameter 5 is equivalent to the Verlet distance 
used in simulations with spherical particles |13j . As shown in the Fig. [2j the 
Verlet method is equivalent to surround the particles by a skin, so that the 
neighbors list consists of all particles pairs whose skins overlap. 

A link cell algorithm [13] is used to allow rapid calculation of this neighbor list: 
First, the space occupied by the particles is divided in cells of side D+S, where 
D is the maximal diameter of the particles. Then the link cell list is defined 
as the list of particles hosted in each cell. Finally, the candidates of neighbors 
for each particle are searched only in the cell occupied by this particle, and 
its eight neighbor cells. 

The neighbor list is calculated at the beginning of the simulation, and it is 
updated when the following condition is satisfied: 

max {Axi + RiAOA > 5. (7) 

l<i<N l J w 

Axi and A6i are the maximal displacement and rotation of the particle after 
the last neighbor list update. Ri the maximal distance from the points on the 
particle to its center of mass. After each update Ax« and A6i are set to zero. 
The update condition is checked in each time step. Increasing the value of 5 
makes updating of the list less frequent, but increases its size, and hence the 
memory used in the simulation. Therefore, the parameter 5 must be chosen 
by making a compromise between the storage (size of the neighbor list) and 
the computer time (frequency of list updates). 
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Fig. 2. Method for identification of neighbor list: the space domain is divided by 
square cells. Then the potential neighbor of the particles are those hosted in the 
same cells, or in the adjacent cells. Each particle has as skin of thickness 5. If the 
skins of two potential neighbors overlap, they are included in the neighbor list. 

Neighbor list reduces the amount of calculations to 0(NM 2 ). Therefore the 
simulations are still very expensive when particles consist of a large number 
of vertices. Further reduction of the number of calculations between neighbor 
particles can be done by identifying which part of a particle is neighbor to 
the other. This idea is implemented as follows: for each element of the neigh- 
bor list, we create a contact list, which consists of those vertex-edge pairs 
whose distance between them is less than ri + r 3 - + 25, where r< and rj are the 
sphero-radii. In each time step, only these vertex-edge pairs are involved in 
the contact force calculations. Overall, neighborhood identification requires a 
neighbor list with all pair of neighbor particles, and one contact list for each 
pair of neighbors. These lists require little memory storage, and they reduce 
the amount of calculations of contact forces to O(N), which is of the same 
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order as in simulations with spherical particles [T3] . 



2.5 Time Integration 



The equations of motion of the system are numerically solved using a four order 
predictor-corrector algorithm [13] . A pseudocode with the basic procedures in 
each time step is shown in Algorithm [TJ The predicted method calculates the 
position (center of mass and orientation) of each particle and its derivatives 
using a Taylor expansion. Next the vertices of the all polygons are updated 
according to the predicted positions of the particles. If the neighbor update 
condition of Eq. [7] is satisfied, the link cell is calculated, and then it is used 
to update the neighbor list and the contact list of each one of its elements. 
Then the contact forces and torques are calculated. Finally, forces and torques 
are used to correct the position of the particles and their derivatives. The 
algorithm is basically the same as the used in polygons, except that the force 
is calculated using Eqs. O Note that the efficiency of the code is based in the 
simplicity of the contact force calculation, and in the fact that the Minkowski 
sum does not need to be calculated during the time integration. 

The parameters of the simulations are a constant stiffness k = 10 7 N/m, grav- 
ity 9 = 10m/s 2 , density a = lkg/m 2 , time step At = 10~ 5 s and Verlet 
distance 5 = Ira. 

Input: state of the particles at time t 
Output: state of the particles at time t + At 
predict position of the particles and its derivatives; 
if neighbor update condition is satisfied then 

calculate link cell; 

update neighbor list; 

update contact lists; 
end 

update vertices of the particles; 

calculate contact forces between neighbor particles; 

apply contact forces to the particles; 

apply gravity forces to the particles; 

correct positions and its derivatives using forces and torques; 

Algorithm 1: One time step of the time integration scheme 
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3 Non-dissipative granular dynamics simulations 



As a first step we present here simulations results of many body conserva- 
tive systems. We investigate the evolutions towards the statistical mechanical 
equilibrium of this simple system. Generalization to dissipative system driven 
by external forces will be presented in forthcoming papers. 

3.1 Energy balance 

The vertex-edge interaction between the particles is given by 

F(V, E) = k8(V, E)N (8) 



where The material parameter k is the stiffness constant, S(V,E) is given by 
Eq. [2] and A" is the unit normal vector: 

N = X-t- (9) 
\\Y-X\\ 



Here X is the position of the vertex V and Y is its closest point on the edge 
E. 

The question which now arises is whether the vertex-edge interaction in Eq. [8] 
leads to a conservative system. Let us multiply the first equation in (jBJ) by r 
and the second one by (p. Next we sum both equations and then sum over all 
particles. After some algebra we get the energy conservation equation: 

Et = E Etj + E ( W + \l^f) = cte. (10) 

ij i 

The first term of this equation corresponds to the potential elastic energy at 
the contacts: 

% = l k (H + E W,Ei) ). (ii) 



The other terms of Eq. (|T0|) are the linear and rotational kinetic energy of the 
particles. We emphasized that the elastic force in Eq. ([3]) belongs to the poten- 
tial energy defined by Eq. ([IT]) , which proves that our model is conservative. 
The simplicity of this force contrasts to the Poschel's model for interacting 



9 



(b) 



(d) 



Fig. 3. Systems obtained from Minkowski sum approach: (a) disks (Vertex + disk); 
(b) rice (segment + disk); (c) peanuts (Polyline + disk) and (d) pebbles (triangle 
+disk) 



triangles where the forces and torques associated to his potential energy 
lead to a much more expensive calculation. 

Other important aspect of this dynamics simulation is the accuracy of the 
numerical solution. The numerical error in the energy calculation is evaluated 
by performing a series of simulations with many non-spherical particles inter- 
action via the elastic force given by Eq. [SJ Each test consists of 400 particles 
confined by four fixed rectangular walls. Each particle occupies an area of 
lcm 2 and the confining area is 46cm x 46.25cm. These dimensions lead to 
a volume fraction of $ = 0.186. Each sample consists of identical particles 
with a specific spheropolygonal shape as shown Fig. [3j disks (point+disks) 
rice (line+disk), peanuts (polyline+disk), and pebbles (triangle+disk). 
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Fig. 4. Time evolution of the total kinetic and potential energy in the non-dissipa- 
tive system. As is expected from the energy conservation, the total energy keeps 
constants, except numerical error wicho are lower than 0.01% thoughout all the 
simulations. 

Initially, each particle has zero angular velocity and a linear velocity of lcm/ s 
with random orientation. Due to collisions, the linear momentum of each par- 
ticle changes and part of it is transferred to angular momentum. FigJH show 
the potential (a) kinetic (b) and total (c) energy of the system. Elastic energy 
has a negligible contribution to the energy budget, as it differs from zero only 
for short times during collisions. Energy conservation is numerically verified 
within a porcentual error of 0.01%. The energy fluctuations are produced by 
time discretization We shall note alse that energy have a trend to grow slowly 
in all samples. 



3.2 Statistical equilibrium 



Simulations with a large number of particles show that elastic interactions 
allow the particles to exchange energy and momentum, whereas its contribu- 
tion to the energy budget is negligible. This property leads us to investigate 
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the existence of a statistical equilibrium in a gas of non-spherical particles. 
It is expected that the system will reach the statistical equilibrium, which is 
characterized by an energy equipartition and a Maxwell-Boltzmann statistics 
for energy distribution [T7] : 

The energy of the system consist of rotational and translational kinetic energy. 
At the beginning of the simulations the rotational kinetic energy is zero, and 
it increases during the simulations due to collision, see left part of the F 
ig. ??. For all the samples, we observe the same stationary regime, where 
the average of rotational kinetic energy reaches the limit of 1/2 of the linear 
kinetic energy. This is in agreement with the theorem of equipartition of energy 
|17j . which states that each quadratic term in the energy should contribute 
the same weight in the mean energy of the system. . It is expected that the 
system reach the statistical equilibrium, which is characterized by a Maxwell- 
Boltzmann statistics for energy distribution [T7j : 

p{E k )dE k = 2^E k /n(kT) 3 exp(-E k /k B T)dE k , (12) 

here E k = ^(mvl + mVy + Iuo 2 ) is the kinetic energy of the particle; T the tem- 
perature; and k B the Boltzmann constant. The mean energy of the particles 
leads to E k = J °° p(E k )E k dE k = \k B T. 

We also calculate the energy distribution of the particles in the stationary 
regime. In order to calculate the energy distribution, we take snapshots of the 
simulations between t = Is and t = 8s distanced by 0.01s. In each snapshot 
the kinetic energy of the individual particles is measured. The histogram of 
the variable e = E k /E k is constructed using 100 identical bins between zero 
and the maximal value. According to Eq. ([121 . the theoretical distribution 
of e must satisfy p(e)de = 2^ej3 3 exp(—/3e)de, where (3 = 3/2. An excellent 
agreement between this theoretical distribution and the simulations data is 
shown in Fig. El Simulations with different non-spherical shapes, which will be 
presented elsewhere [H], show that energy distribution for all samples collapse 
onto the theoretical expected value. It is also shown that the relaxation time 
for the statistical equilibrium is very sensitive to the degree of non-sphericity 
of the particle. 



4 Dissipative granular dynamics simulation 

Here we present biaxial tests simultions with circular and non-circular parti- 
cles. Usually, the granular assemblies are compacted and loaded within a set 
of confining walls. These walls act as boundary conditions, and can be moved 
by specifying their velocity or the force applied on them. The response of the 
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1 2 3 4 5 

time (t/tc) x10 4 

Fig. 5. Time evolution of the total linear and rotational kinetic energy of the parti- 
cles. Eq is the initial value of the total kinetic energy. The horizontal lines correspond 
to the expected value by the equilibrium statistical mechanics. 

walls can be used to calculate the global stress and strain of the assembly. 

The interaction of the spheropolygonspolygons with the walls is modeled here 
by using a simple visco-elastic force. First, we allow the polygons to penetrate 
the walls. Then, for each vertex of the polygon a inside of the walls we include 
a force 

Confining walls can be used to generate samples with different void ratios. 
Starting from a very loose packing, the sample is compacted by applying 
a centripetal gravitational field to the particles and on the walls, oriented 
to the center of mass of the assembly. Then the sample is subjected to an 
isotropic compression until the desired confining pressure is reached. In order 
to generate dense samples, the interparticle friction is set to zero during the 
construction. The loose samples are created taking damping coefficients 100 
times greater than those used in the test stage. Samples with void ratio ranged 
from 0.128 to 0.271 can be achieved with this method |19j . 

We have investigated shear deformation of granular samples with different 
initial void ratios [19]. Shear bands are observed in dense samples, whereas 
they seems to be absent in loose ones, they share some common properties of 
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v peanuts 




Fig. 6. energy distribution N(e) for the particles, where e = Ek/E^-The line corre- 
sponds to the best fit n(e) = 2/3 ef3/it exp(— /3e), with = 1.5. 



the shear bands in real granular materials, such as their characteristic reflec- 
tion when they reach the boundary wall. Shear band orientation lies between 
the Roscoe angle and Mohr-Coulomb solution, as in most of the experimental 
data. 



For large shear deformations all samples reaches the critical state, which is 
independent on their initial density. Once the samples reaches this state, they 
deform at constant void ratio and coordination number [19] . The evolution of 
the deviatoric stress exhibits fluctuations around the residual strength. Abrupt 
reduction of the stress results from the collapse of force chains, as shown the 
Fig.??, collapse of force chains makes the sample to approach and retreats 
unstable stages. A similar behavior is observed in glass bead samples [20] and 
packings of glass spheres |21j . Experimental biaxial tests show evidence of 
dynamic instabilities at the critical state [22] . Erratic slip-stick motion at the 
critical state is interesting, owing to its potential analogy with earthquake 
dynamics [2"5] . 
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Fig. 7. Granular packing obtained with disks, rice, peanuts and rice 
5 Performance 



Lastly, We compare the efficiency of the many-body simulations of systems 
consisting on disks, spheropolygons and clums of spheres. Each spheropoly- 
gons correspond to a particle with complex shape, and it consist on 62 ver- 
tices. The clump of spheres represents the same complex shaped particles, and 
it need 726 particles. The macroparticles they are simulating by summing 
up the contact forces between the constituting disks, and updating all the 
disks of each particle according with its current position. The performance 
of the simulations is estimated by running different processes in a Pentium 4, 
3.0GHz, and calculating the Cundall number is each one of them. This number 
is the amount of particle time steps executed by the processor in one second, 
which is calculated as c = N T N/T C pu, where N T is the number of time steps, 
N is the number of particles and Tqpu is the CPU time of the simulation. 
Fig. shows the Cundall number versus the number of particles for the three 
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Fig. 8. Cundall number versus the number of particles, in simulations with disks, 
spheropolygons and clumpy particles. 

cases. When the number of particles is between N = 100 and N = 1000 the 
Cundall number is approximately constant. This constant is around 100, 000 
for disks, 2, 000 for spheropolygons, and 50 for clumpy particles. Therefore 
the simulations with spheropolygons, although slower than simulations with 
disks, are much faster than simulations with clumpy particles. This is because 
each time step needs to update the position of 62 vertices of the shperopoly- 
gons whereas it needs to update the position of the 726 disks in the case of 
the clumpy particles. Therefore simulations with spheropolygons are more ef- 
ficient than those ones with clumps of disks, because the former ones require 
less elements to represent the particle shape. 



6 Concluding remarks 



The method presented here provide an energy balance equations and a wide 
range of particle shape representations, including non-convex particles and 
tunable grain roundness. This paper shows that modeling interacting parti- 
cles using spheropolygons has several advantages with respect to other existing 
particle-based models: i) The possibility to model non-convex particles; ii) a 
realistic representation of the surface curvature of particles; iii) guaranteed 
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(a) 
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strain (A H/Ho) 



Fig. 9. Stress versus strain (a) and void ration versus strain (b) for different particles 
shape, in biaxial test simulations. The dotted line in (a) represents p/k n , where p is 
the applied pressure on the latteral walls. 
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compliance with physical and statistical mechanical laws; iv) balance between 
accuracy and efficiency Benchmark tests prove energy conservation with an 
error below 0.0001%. Simulations with many particles verify the Maxwell- 
Boltzmann distribution and the principle of energy equipartition. The com- 
putational efficiency is compared to simulations with disks and clumps of 
disks. Simulations with disks are around 50 times faster than simulation with 
spheropolygons. However, the speed of the simulation is 40 times faster that 
simulations with clumpy particles. 

THe method overcome also two main difficulties existing in previous DEM 
developments. 

(1) The interaction between polygons using the overlapping area is difficult 
to generalize in 3D, because the overlap between polyhedrons is much more 
difficult to evaluate. 

(2) The elastic force used in this work does not belong from a potential, 
so that this model does not provide an equation for energy balance. In the 
investigation of fault zones, the energy balance is required to determine the 
energy budget in earthquakes. 

The model is still very simple, but extensions to more complex interactions 
and 3D simulations are achievable in the near future. Cohesive and frictional 
forces can be incorporated by introducing internal variables in each vertex- 
edge contact. These variables account for elastic deformation at the contacts 
and they must be updated in each time step according to the sliding condi- 
tions or breaking criterion. 3D modeling using spheropolyhedra requires elastic 
forces similar to Eq. ([3]), where the sum is over all vertex-face and edge-edge 
interactions. Special attention is required for the case of two parallel edges in 
contact. This case lead to a non- uniqueness in the selection of contact points. 
This need to be resolved to get a physical correctness in the torque calculation. 
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